Symmetry-projected variational approach for ground and excited states of the 

two-dimensional Hubbard model 
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We present a symmetry-projected configuration mixing scheme to describe ground and excited 
states, with well defined quantum numbers, of the two-dimensional Hubbard model with nearest- 
neighbor hopping and periodic boundary conditions. Results for the half-filled 2 x 4, 4 x 4, and 6x6 
lattices, as well as doped 4x4 systems, compare well with available results, both exact and from 
other state-of-the-art approximations. We report spectral functions and density of states obtained 
from a well-controlled ansatz for the (N £ ± l)-electron system. Symmetry projected methods have 
been widely used for the many-body nuclear physics problem but have received little attention in 
the solid state community. Given their relatively low (mean-field) computational cost and the high 
quality of results here reported, we believe that they deserve further scrutiny. 

PACS numbers: 71.10Fd, 21.60.-n 
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I. INTRODUCTION 

Since the discovery of high-Tc superconductivity^ 
there has been a growing interest in the properties of cor- 
related two-dimensional (2D) electronic systems. 2 Within 
this context, the Hubbard model 3 - has received a lot of 
attention since it is considered one of the simplest mod- 
els still containing the relevant physics^ Renewed inter- 
est in the Hubbard Hamiltonian also comes from recent 
experiments 5 -^ with cold fermionic atoms in optical lat- 
tices which open the possibility for direct simulations of 
the model with lattice emulators. 7 Hubbard-like models 
are also relevant to describe electronic properties within 
the active research field of graphene£ 

The repulsive Hubbard Hamiltonian is a very inter- 
esting model in theoretical physics. On the one hand, 
neither its hopping (one-body) nor its on-site interaction 
(two-body) terms favor any interesting magnetic order- 
ing. On the other hand, when both of them combine into 
the full Hamiltonian a rich variety of interesting phenom- 
ena is displayed, for example, correlation-driven metal- 
insulator transitions^ ferromagnetism,— deviations from 
the standard Fermi-liquid results*** long-wavelength col- 
lective modes** and spatially inhomogeneous phases^- 3 - 
The dimensionality of the model also challenges the the- 
oretical tools at our disposal. Exact analytical solu- 
tions exist in the one-dimensional (ID) case** whereas 
the present knowledge of the basic quantum mechani- 
cal properties of the 2D Hubbard Hamiltonian relies, to 
a large extent, on numerical techniques applied to the 
Hamiltonian itself or to its strong coupling approxima- 
tions, i.e., the t-J, t-J* and Heisenberg models^ ' 15 ' 16 In 
particular, for the case of the full 2D Hubbard Hamilto- 
nian, a very efficient Lanczos algorithm*** based on the 
classification of all the irreducible representations of the 
space group, has allowed systematic studies in the 4x4 
lattice. 



Going beyond the present limits of exact diagonal- 
ization (ED) techniques requires a truncation strategy. 
A key issue is then how to truncate the model space 
while still being able to retain the most important de- 
grees of freedom relevant for the description of a partic- 
ular ground and/or excited state. Nowadays there are 
several methods at our disposal, some of them already 
heavily used to study ID and 2D Hubbard models with 
variable degree of success. One that has been used with 
great success is the Quantum Monte Carlo** - — (QMC) 
approach. Another is the density matrix renormalization 
group** - — (DMRG) scheme that represents a very pow- 
erful and general decimation prescription. Currently, the 
DMRG algorithm is understood as an energy minimiza- 
tion within a class of low entanglement wavefunctions 
known as matrix product state s 24 ' 25 (MPS) establishing 
an exciting link with quantum information perpectives^SS 
A very flexible entanglement encoding is also provided 
by the rapidly expanding research area of tensor network 
states*^.** (TNS). 

Variational principles also offer very powerful meth- 
ods to study Hubbard-like models. For example, the dy- 
namical variational principle ) 30 ' 31 expressed in the lan- 
guage of Green's functions and self-energies^* provides 
us with the variational cluster approximation 33 (VCA), 
the dynamical impurity approximation 34 (DIA) and the 
dynamical mean field theory** (DMFT). Within this 
context, the self-energy-functional theory 3 - (SFT) has 
emerged as a conceptual framework in which the VCA, 
DIA and DMFT, as well as several extensions of them, 
can be specified by the choice of a reference system. In 
particular, the cluster extensions to DMFT have pro- 
vided important insights into the physics of the 2D Hub- 
bard model in aspects such as the Mott-Hubbard transi- 
tion, the pseudogap in doped systems, and the phase dia- 
gram itself ,2L2& DMFT and its cluster extensions are par- 
ticularly valuable as they have been shown to be comple- 
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mentary to finite size simulations ; i 39 — including ours. 
Here, we also refer the interested reader to recent work 42 
where a hierarchy of truncated configuration interaction 
(CI) expansions has been considered as a solver for quan- 
tum impurity models and DMFT. 

In the present work, we explore an alternative av- 
enue not only to describe ground state properties of the 
2D Hubbard model but also to access excitation spectra 
which represent a basic fingerprint of quantum mechan- 
ical correlations in the considered lattices. A first step 
in this direction, based on symmetry-projected configu- 
ration mixing ideas originally employed in microscopic 
nuclear structure theory^ 3 - was undertaken for the ID 
Hubbard model^i and is extended in the present work 
to the 2D Hubbard Hamiltonian with periodic boundary 
conditions (PBC). 

For a given single-electron space, we construct the most 
general unitary Hartree-Fock (HF) transformation^ 4 ^ 
Since this HF-transformation mixes all the spin and lin- 
ear momentum quantum numbers of the single-electron 
basis states, the corresponding Slater determinant delib- 
erately breaks the original spin and translational sym- 
metries of the 2D Hubbard Hamiltonian. Therefore, as 
such, our symmetry-broken Slater determinant can be 
considered as a convenient mean-field starting point en- 
larging the space of trial wave functions We restore 
the broken translational and spin symmetries with the 
help of linear and angular momentum projection oper- 
ators. This symmetry restoration recovers the multi- 
determinantal character in our trial state keeping good 
spin and linear momentum quantum numbers. The Ritz 
variational principl o 45 ' 46 is then applied to the projected 
energy, i.e., ours is a variation-after-projection (VAP) 
scheme. This procedure provides us with the optimal 
(variational) representation of a ground state, with well 
defined spin and linear momentum quantum numbers, 
via a single symmetry-projected configuration. Our VAP 
scheme is also very close in spirit to Projected Quasiparti- 
cle Theor y 47 ' 48 (PQT) and is related to other variational 
approaches , 49 i 50 

In order to describe excited states with well defined 
quantum numbers, we construct a truncated basis con- 
sisting of a few (orthonormalized) symmetry-projected 
states throughout a chain of VAP calculations. This can 
be easily done, still with low computational cost, due to 
the simple structure of our projected wave functions. Fi- 
nally, a further diagonalization of the 2D Hubbard Hamil- 
tonian is performed within such a basis. With this config- 
uration mixing procedure we may account, in a similar 
fashion, for additional correlations in both ground and 
excited states. In addition, our theoretical framework can 
be used to study important dynamical properties of the 
2D Hubbard Hamiltonian like spectral functions ^ 15 ' 32 

In this paper we have three main goals. First, we 
present the methodology of a VAP configuration mix- 
ing scheme, originally devised for the nuclear many-body 
problem, but not yet explored for the 2D Hubbard model. 
Therefore, in Sec. |TT]we introduce our theoretical formal- 



ism. Symmetry restoration is described in Sec. Ill Al while 
our configuration mixing scheme is outlined in Sec. IHBI 
For the reader's convenience, the key ingredients of our 
approximations are stressed in these two sections while, 
to make our presentation self-contained, more technical 
details can be found in appendices IA1 and iBl respectively. 
Our second goal is to show how our theoretical frame- 
work can be used to access the spectral weight of states 
with different linear momentum quantum numbers. To 
this end, the computation of hole and particle spectral 
functions is briefly described in Sec. Ill CI and more de- 
tails are given in appendix[Cj Our third goal is to test the 
performance of our approximation for a selected set of il- 
lustrative examples. The results of our calculations for 
the half-filled 2 x 4, 4 x 4 and 6x6 lattices are discussed 
in Sec. Mil There, we pay attention to the properties 
of ground and excited states but also discuss hole and 
particle spectral functions as well as the corresponding 
density of states (DOS). In addition, in the case of the 
4x4 lattice, we consider doped systems with 14 and 15 
electrons. Finally, Sec. II VI is devoted to the concluding 
remarks and work perspectives. 

II. THEORETICAL FRAMEWORK 

In what follows, we describe the theoretical framework 
used in the present study. First, symmetry restoration 
and configuration mixing are presented in Sees. Ill Al and 
IIIBI The computation of spectral functions is briefly de- 
scribed in Sec. Ill CI 

A. Symmetry restoration for the 2D Hubbard 
model 

We consider the following one-band version of the 2D 
Hubbard Hamiltonian 3 

H H ub = -tJ2 ( £ j+x<r e ^ + 4+y<? £ i° + h - c ) 

+ u T, d lAiwn (!) 

where the first term represents the nearest-neighbor hop- 
ping (t > 0), with unit hopping vectors x = (1,0) and 
y = (0, 1), and the second is the repulsive on-site inter- 
action (U > 0). The operators cj a and c^ a create and 
destroy a particle with spin-projection a = ±1/2 (also 
denoted as a =^, 4,) along an arbitrary chosen quantiza- 
tion axis on a lattice site j= (j^ , jy) . They satisfy the 
usual anticommutation relations for fermion operators*^, 
Here, and in what follows, the lattice indices run as j x — 
1,...,N X and j y = 1, . . . , N y with N x and N y being the 
number of sites along the x and y directions, respectively. 
The total number of sites is given by N S i tes — N x x N y . 
We assume PBC, i.e., the sites Ni + 1 and 1, with i=x,y, 
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FIG. 1: (Color online) The energy spectrum, obtained via Eq. (|21[) . for the half-filled 2x4 lattice at U=4t is shown in panel 

a) . This spectrum can be hardly distinguished from the one obtained using an exact diagonalization (ED). Therefore, in panel 

b) the absolute errors are plotted for each of the predicted 120 solutions. For more details, see the main text. 



are identical. Furthermore, we assume a lattice spacing 
A=l. 

Next, we apply the 2D Fourier transform 

to obtain operators with momentum k Q = (fc Q;c ,fc Qy ) = 

(^T> ^wf) ■ Thc Hamiltonian (JTJ) can be easily written 
in terms of these new operators. The quantum numbers 
Oj, with i=x,y, take the allowed values 



P =VD* c 1 (4) 

u a / j " aa, a'-aa V*) 

tx<7 

where I? is a general 2N S i tes x 2N S i tes unitary 
transformation ! 45 ' 46 In Eq. ^ a is a shorthand notation 
for the set (a x ,a y ,cr a ). The transformation Q mixes all 
the linear momentum quantum numbers as well as the 
spin projection of the states ©. As a consequence, \T>) 
deliberately breaks rotational (in spin space) and transna- 
tional invariances. To restore the spin quantum numbers 
we explicitly use the projection operator 



inside the Brillouin zone (BZ)^i Equivalently, they can 
take all integer values between and Ni — 1. 

In the HF-approximation, the ground state of an N e - 
electron system is represented by a Slater determinant 
\V) = m which the energetically lowest N e 

single-electron states (hole states h, h , . . . ) are occupied 
while the remaining 2N S i te s — N e states (particle states 
p, p , . . . ) are empty. The HF-quasiparticle operators 
are given by 



2^-1-1 f 

pL> = -^J dM> s £,(n)Rs(n) (5) 

where i?s(fi) = e~ %aSz e~ l/3Sv e~ ijSz is the rotation oper- 
ator in spin space, f2 = (a, /3, 7) stands for the set of Eu- 
ler angles and f^s' are Signer functions.— The form 
([5]) has been frequently used for total angular momentum 
projection in nuclear physics i 43 ' 45 This form has also been 
adopted in the study of the ID Hubbard model^ 4 - and 
more recently within PQT in quantum chemistry . 47 ! 48 

The linear momenta, and k^ , are restored with the 
projector 
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FIG. 2: (Color online) The DOS Af(oj) [Eg.lpO))] for the half-filled 2x4 lattice at U=4t is plotted in panel a) as a function of 
the shifted excitation energy u) — U/2 (in t units). Results have been obtained by approximating the (N e ± l)-electron systems 
[Eqs. (|23[) and l|27[l] with pit = 1 (red) and nr = 5 (blue) Slater determinants out of Sec. Ill Al As can be observed from panel 
b) the DOS obtained with exact diagonalization (ED) and the one obtained using tit = 5 HF-transformations can hardly be 
distinguished. The hole (blue) and particle (black) spectral functions, computed with tit = 5 HF-transformations, are plotted 
in panel c). A Lorentzian folding of width F=0.05t has been used. 
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where P 



OLCJ 



(6) 



T c aCT is the generator of 



the considered lattice translations. Note that this oper- 
ator neither has vector character nor corresponds to the 
true linear momentum operator. It is associated with 
the quasi-momentum resulting from translational invari- 
ance of the lattice. We will refer to it, however, as lin- 
ear momentum for simplicity. The projector ^ repre- 
sents the 2D limit of the general operator restoring Galilei 
invariancei 43 ' 53 ' 54 Note that, at variance with atomic nu- 
clei, lattice systems can have solutions with linear mo- 
menta different from zero. 

In what follows, we introduce the shorthand notation 
= (5, £) for the set of (symmetry) quantum numbers 
(S, &,,£„), i.e., P^,C(£) = We then use the fol- 

lowing symmetry-projected wave function 



|P;6;E) = 



s 

E'=-S 



ftp? 



EE' 



\V) 



(7) 



where /®, are variational parameters. Note that, through 

the action of the projection operator P®^ , the multi- 
determinantal character of the state characterized by the 
quantum numbers 6 and S is recovered and written in 
terms of the quantum numbers G and E4£ In practice, 
the integration over the set of Euler angles in Eq.© is 
discretized. For the integrals in a and 7 we have used 
8 grid points whereas for the /3-integration we have used 
16 points. Therefore, the total number of grid points to 
be used for the projection operator P® £ , is 1024 x N sites . 

For a given symmetry O, the energy (independent of 
£) associated with the state (JTJ) 
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(8) 



is given in terms of the (25 + 1) x (25 + 1) Hamiltonian 
H^, = (V\H Hub P^,\V) and norm = (V\P^,\V) 
matrices (see appendix |XJ . It has to be minimized with 
respect to the coefficients / and the HF-transformation 
T>. The variation with respect to the mixing coefficients 
yields the following generalized eigenvalue equation 
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(9) 
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FIG. 3: (Color online) Occupation numbers [Eq. (|26|l ] of the basis states in the ground state of the half-filled 2x4 lattice are 
plotted for various U strengths. 



with the constraint / e tA/" / e = l 2 5 + i ensuring the or- 
thogonality of the solutions. The unrestricted minimiza- 
tion of the energy ([5]) with respect to the underlying 
HF-transformation T> can be carried out via the Thou- 
less theorem i 43 i 45 The corresponding variational equa- 
tions assume the form 



with 



Lr.„ 
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/°t (/C e - E e U e ) f 



(10) 



ph 



(ii) 



B. Symmetry- projected configuration mixing for 
the 2D Hubbard model 



An accurate description of excited states in a many- 
fermion system is much more difficult even when one is 
usually interested in just a small fraction of the low-lying 
spectrum. Here, the main difficulty in the optimization 
of excited states is ensuring orthogonality among them 
and with respect to the ground state. For this, we simply 
use a Gram-Schmidt orthogonalization. Our goal in this 
section is to construct, throughout a chain of VAP cal- 
culations, a basis of a few (orthonormalized) states with 
well defined quantum numbers 0. 

Suppose we have generated a ground state solution 
Icj) 1 ) = |X>;6;£) out of Eqs. © and CEU]) in Sec. UTAH 
Then wc write the first excited state wave function as 



Here, the N e x N e and (2A sites - N e ) x (2N sltes - N e ) 
matrices Lq and Mq are obtained via the Cholesky 
decompositions 1 43 i 53 The particle-hole kernels K^^j 1 = 

{V\H Hub P^,blb h \V) and 1& = (V\P^,blb h \V) are 
given in appendix [B] It should be stressed that, for a 
given symmetry O, we only retain the energetically lowest 
solution of Eqs. © and (TTU1) . Both the HF-transformation 
T> and the mixing coefficients / e are essentially complex, 
therefore one needs to minimize n var = 2(2N S i tes — N e ) x 
N e + AS real variables. We use a quasi-Newton method 
for such a minimization! 55 ' 56 The variational procedure 
already described is known in nuclear structure physics as 
the VAMPIR (i.e., Variation After Mean field Projection 
In Realistic model spaces).— Note, that particle number 
projection^ 5 - is not carried out in the present study since 
the considered Slater determinants conserve the number 
of electrons. 



(12) 



where \c/> 2 ) has a form similar to Eq.([7]) but with differ- 
ent coefficients / 2e and underlying HF-transformation 
T> 2 . The label 2 distinguishes them from the ones (i.e., 
/ ie and V 1 ) corresponding to the reference ground state 
we already have. Both f3f and /3| can be obtained by 
requiring that {<f> l \ip2) — and {^2^2) = 1- They are 
given in terms of the projector (i.e., S\ = Sf) 
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as follows 
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FIG. 4: (Color online) The energy spectrum, obtained via Eq. (|21l) . for the half-filled 4x4 lattice at U=4t is shown in panel a). 
In panel b), the excitation energies from the ground state to the lowest-lying S=l and S=2 states from panel a) are plotted as 
functions of the linear momentum quantum numbers F = (0, 0), Ri = (1, 0), P = (2, 0), R2 = (2, 1), Q — (2, 2), and R3 = (1, 1), 
respectively. In addition to U=4t (blue boxes), results for U=0t (red diamonds) are also included for comparison. 
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The first excited state is obtained varying the energy 
functional for (TT2|) with respect to f and V 2 . For the 
second excited state, we introduce a new state |</> 3 ), again 
with the same form as in Eq.([7]), and write 



= mm + mm + put 3 ) 



(15) 



with coefficients (3f, /?| and /3f such that \tpz) is orthogo- 
nal to the previous solutions \ipi) = [Eq.0] and \y>2) 
[Eq. (TT2")) ] as well as ((p 3 \ip s ) = 1. The second excited state 
is obtained varying the energy functional for (fT5jl with 
respect to / 3e and T> 3 . Let us have a more general situ- 
ation in which, by successive variation, i = 1, . . . , m — 1 
orthonormalized solutions (for example, \ipi) and \ip2)) 



m — l 



(17) 



with \<j} m ) having again the form ([7]). Requiring orthonor- 
malization with respect to all the previous m — l solutions 
(□U) the coefficients /3™ and /3™ in Eq.(H7J) read 



mn = 
r'm 



m — l 

E 



l-S 



m — l 



l<F) 



-1/2 



-B r ' 

n 



in terms of the projector (i.e., S m -i — S^-i) 



E 



The energy for the state ([T7|) takes the form 



(18) 



(19) 



Ei^ 



(16) 



j'=i 



are already at our disposal. Each of the states \<j>>) in (TT6"]) 
has the same form as ([7]) . One then writes the ansatz for 
the mth state wave function (for example, 1 773) ) as 



!!)(-) 



(20) 



with kernels 'H^ 1 ® and Af™® accounting for the fact 
that m-1 linearly independent solutions have been re- 
moved from the variational space. Their expressions 
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FIG. 5: (Color online) The DOS Af(ui) [Ed. 1(50)1] for the half-filled 4x4 lattice at U=4t is plotted in panel a) as a function of 
the shifted excitation energy ui — U/2 (in t units). Results have been obtained by approximating the (N e ± l)-electron systems 
[Eqs. I|23p and (|27[l ] with nr = 5 HF-determinants. Hole (blue) and particle (black) spectral functions, are displayed in panel 
b). A Lorentzian folding of width r=0.2t has been used. 



are slightly more involved^ than the ones required in 
Eq.© but still straightforward. They require the knowl- 
edge of the symmetry-projected matrix elements between 
two different Slater determinants \T> 1 ) and \T> k ) (see ap- 
pendix The variation of the energy (|20p with respect 
to f m& yields an equation similar to (JSJ) with the con- 
straint /m©tjV" me / me = 12S+1- The unrestricted mini- 
mization of the energy ([20]) with respect to V m , via the 
Thouless theorem, leads to variational equations similar 



to mi but with kernels K%Z <pn and K" l °; ph that re- 
quire symmetry-projected particle-hole matrix elements 
between two different Slater determinants |2? 4 ) and \D k ) 
(see appendix IB1 . 

The procedure outlined in this section is known in nu- 
clear structure physics as EXCITED VAMPIR. 43 It pro- 
vides a (truncated) basis of m (orthonormalized) states 
with a well denned symmetry O, still keeping low 
computational cost. This is doable due to the simple 
structure of the projected states defining such a basis in 
combination with a fast minimization schcmei 55 ' 56 Our 
method can also be extended to use general Hartrcc- 
Fock-Bogoliubov (HFB) transformations ! 43 ' 45 ' 47 How- 
ever, this requires an additional projection of the particle 
number, which increases the numerical effort by about 
one order of magnitude and has hence not been used in 
the present paper. 

It should be noticed that the ground state \tpi) [Eq.©] 
is written as a projection operator acting on a single de- 
terminant, the first excited state | ^2) [Eq. (TT2")) ] as a pro- 



jection operator acting on two determinants, and so on. 
Because this allows excited state wave functions to be 
described at a higher level of quality than is the ground 
state wave function, our final step is to diagonalize the 
2D Hubbard Hamiltonian in the basis of the states \<fij). 



(21) 



For ground and excited states, the resulting wave func- 
tions 



(22) 



may account for more correlations than the description 
based on a single symmetry-projected configuration dis- 
cussed in Sec. Ill Al In the present work, as a first step, 
we have restricted ourselves to test the performance of 
our approximation with m=5 (orthonormalized) states. 
As we will see, this turns out to be a reasonable start- 
ing point for, at least, a qualitative description of the 
considered lattices. 

An interesting issue is the evolution of the energy of 
each state with the number m of transformations in- 
cluded in the prescription described in this section. We 
observe that, for the lattices considered in the present 
study, the energy of the ground and the first couple of 
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excited states remains unchanged when m goes from 1 
to 5 (the changes in the energy per site are of the or- 
der 10 ). This is partly because the main correlations 
have already been accounted for with a single symmetry- 
projected determinant. Therefore, the excited configu- 
rations obtained constitute reasonably good approxima- 
tions to the true excited states of the considered system. 
We produce m=5 symmetry-projected determinants in 
order to obtain the low-lying spectrum. For the sys- 
tems considered in this work, these states turn out to 
be weakly coupled through the Hamiltonian. However, 
this cannot be anticipated a priori and the diagonaliza- 
tion Eq. (|2~lj) should always be carried out. Preliminary 
results for larger square lattices (i.e., 8x8 and 10 x 10) 
as well as for other Hamiltonians (i.e., the t — t — U 
and t — t —t—U Hubbard models) indicate that there 
are cases in which the diagonalization Eq. ([2Tj) brings a 
sizeable amount of additional correlations. 



C. Hole and particle spectral functions 

Let us assume that for an even number N e of elec- 
trons we already have the ground state wave function 
IX? 1 ; ©° ; = 0), out of the calculations described in 
Sec. Ill Al Since for all the considered lattices with 
an even number of electrons the ground state has spin 
S=0, but not neccessarily linear momenta zero, we write 
its quantum numbers as 9° = (0, £ ). Usually, spec- 
tral functions are computed within a Green's function 
perspective^ The key point is then to approximate the 
ground states of the (N e ± l)-electron systems by a suit- 
able ansatz. In the present study, we approximate^ the 
ground state of the (-/V e -l)-electron system, with the sym- 
metr y e- = (5 = l/2,r), by 



= E f? h 



„ hl P^b h (V*W) (23) 



iha 



where the index i runs as i — l,...,nr, the hole in- 
dex h as h = 1, . . . , N e and a = ±1/2. In Eq.(|2"3"]l. we 
write bh{'D t ) to explicitly indicate that holes are made 
on tit different Slater determinants. The determinants 
II? 1 ) and \T> Z ) correspond to the ground and lowest en- 
ergy (i = 2, . . . , nr) states obtained for the iV e -electron 
system out of the calculations described in Sec. Ill Al In 
the present study, we have restricted ourselves to a max- 
imum of tit — 5 HF-transformations. The coefficients 
/ e in Eq. ([23l) are obtained by solving the equation 



(24) 



(H e ~ - Af e ~) f e ~ =0 



that yields 2riTN e hole solutions hi with ener- 
gies E^ . With all the previous ingredients, 
one can compute^— the hole spectral function as 
S hl (C,Se hl ) = |<fc i; e-||^o_ € - HI? 1 ; 0°>| a in terms of 
the reduced matrix element 



ir^e ) = -- 



i 



8TT 2 N sites Y (E>i|pe°|X>i 

E &E e ^ j / d0i) -'* (n)( - 1)1/2 " a ' x 



ihh aa 



f) 1 * 



(25) 



where = (k^-,k^-J = ( 2 ^ x , y j . The indices 
i, h, h and a, a run as in (|23[) . £~ and £~ run as in 
Eq.© and 5e hl = E e ° ~ E® . Details for the computa- 
tion of the kernels H and A/" e in Eq.(|24|) as well as 

r i _1 

Xj} h (il,j) and n ll {Q, j) in Eq.([25]) can be found in 

appendices [O and [Al respectively. The occupation num- 
ber n(£ _ ) of a basis state @ in the ./V e -electron ground 
state can be computed as 



2n T N c 
hi = l 



(26) 



The (iV e +l)-elcctron system, with the symmetry 8 + = 
(S = 1/2, £ + ), is approximated by^ 

|p i; e+;a) = Y^C'^K^W) (27) 

ipa 

where the index i runs again as in (|23|) . The particle index 
p takes the values p = N e + 1, . . . , 2N S i tes and a — ±1/2. 
In this case, the coefficients g e are obtained by solving 
the equation 



(28) 



that yields 2tlt{2N sites — N e ) particle solutions p\ with 
energies E® . The particle spectral function is then 
written as S Pl (£ + , 5e pi ) = |(pi; e+||c € +_ $ o||X> 1 ; 9°)| 2 in 
terms of the reduced matrix element 



<Pi;e+ 



i| P_e°l2?i\ 



/ <P \* oo 

E ^ E ^ J xutJfvK, (n,i)x 



ipp a a 3 



(29) 



where, in this case, k^ = (jt^+,k^+^j = ( , \ . 
The indices i, p,p and a, a run as in (|27jl. and 
run as in Eq.Q and 5e Pl — E®^ — E e ° . Details for 
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FIG. 6: (Color online) The same as Fig[S]but for U=8t. The shapes of the DOS as well as the spectral functions for momenta 
{it, 0), (7r/2,0) and (0,0) are qualitatively similar to the ones obtained using Lanczos calculations.— 



the computation of the kernels "H 0+ and A/" e+ as well as 



,{ft,j) in Eq. (j2T))) can be found in appendix [Cl 
Finally, the DOS can be computed as 



(30) 



where the indices hi and pi are absorbed into the contin- 
uous variable lu. Due to the finite size of the system the 
spectral functions consist of a finite number of S func- 
tions with different weights. Therefore, we introduce an 
artificial width T for each state using a Lorentzian. In all 
cases our DOS is normalized to 2 x N sites . 



III. DISCUSSION OF RESULTS 

In this section, we discuss the results of our study. We 
have considered the 2x4 half-filled lattice as a prototyp- 
ical system where one can obtain the full spectrum by 
means of ED. This allows us to callibrate our approx- 
imation not only for ground state properties but also 
for excited states. Next, we have considered the well- 
studied half-filled 4x4 lattice, which constitutes the 
largest square lattice for which exact ground state ener- 
gies are available in the literature. Other approximation 
schemes have also been tested for this lattice in previous 
works. Results have already been published for doped 
systems with 14 and 15 electrons in this lattice, which 
motivated us to also perform calculations for them in the 



present study. Last, we consider the half-filled 6x6 lat- 
tice as a prototype of a system where ED is no longer 
feasible. Many of the results to be discussed in what fol- 
lows correspond to U=4t taken as a representative on-site 
repulsion for which studies are available. Nevertheless, 
let us stress that our formalism can be used for any 2D 
Hubbard hamiltonian of the form ([1} with arbitrary U 
and/or t values. 



A. The square 2x4 lattice 

Let us start by considering the rectangular 2x4 lat- 
tice. The first five solutions obtained at half-filling via 
Eq. (l21[) . for each of the linear momentum quantum num- 
bers (0,0), (0,1), (0,2), (0,3), (1,0), (1,1), (1,2) and (1,3) 
and the spins S=0,1, and 2, are plotted in panel a) of 
Figfflfor U=4t. The first excited state corresponds to a 
O = (1,1,2) configuration [with linear momenta (tt, it)]. 
The energies e® of the 120 solutions shown in the figure, 
have been compared to the ones obtained using an ED. 58 
The comparison reveals that both spectra follow the same 
qualitative trend and can hardly be distinguished. There- 
fore in panel b) of the same figure, we have plotted the ab- 
solute errors e a bsoi. — E exact — E for each of the predicted 
120 states. Our approximation fairly reproduces the ex- 
act ground state energy -10.2529t for this system. For all 
the 40 S=0 and S=l solutions considered the absolute er- 
rors remain very small, the largest deviation being 0.047t 
for the second state with the symmetry O = (1,0,0). 
The previous results are encouraging if one takes into 
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FIG. 7: (Color online) Energy spectrum, obtained via Eq.lJ^TJ), for the 4x4 lattice with N e = 15 electrons at U=4t. 



account that, even for this relatively small lattice, the 
number of variational parameters in our approximation 
n var (S = 0,£ x ,€ v ) = 128 and n var (S = l,€ x ,€ v ) = 132 
is about half of the dimensions (S 1 = 0, £ X) £ y ) — 221 
and nRH(S — l,6r,f y ) — 294 of the restricted Hilbert 
spaces. On the other hand, n var (S = 2,£ x ,£ y ) = 136 is 
larger than tirh{S — 2,£ x ,£ v ) = 90 and therefore our 
solutions reproduce the ED ones for S=2 states. 

In panel a) of Fig. [2j we have plotted the DOS M{oj) 
[Eq.dSD])] for the half-filled 2x4 lattice at U=4t. The cal- 
culations have been carried out by approximating the (N e 
± l)-electron systems [see Eqs. (|23]l and (f2T]> ] with nr=l 
(red curve) and n^—h (blue curve) HF-transformations 
along the lines described in Sec. Ill CI We have intro- 
duced a shift equal to the chemical potential at half-filling 
(/Jo = U/2) so that the DOS in FigJ2] appears to be sym- 
metric around w-U/2=0. This convention, i.e., to plot 
DOS and spectral functions vs. w-U/2 will be adopted in 
the rest of the paper. The DOS shows the Hubbard gap, 
Ah = U/2 = 2t, characteristic of finite size lattices. We 
note, however, that previous studies within the frame- 
work of the dynamical cluster approximation (DCA) have 
shown that the gap is preserved at sufficiently low tem- 
peratures even in the thermodynamic limit (TDL)^2r— 
On the other hand, the nonperturbative study of Ref. 59 
has concluded that for the half-filled Hubbard model the 
gap persists for any finite value of the on-site repulsion 
U, the only singular point being U=0t. 

From panel a) of FigJUone realizes that, even for this 



small lattice, the fine details of the energy distribution of 
J\f(uj) can only be obtained using a larger number tit— 5 
of HF-transformations to describe the (N e ± l)-electron 
systems. Using tit— 5 transformations, Eqs. dMl) and ([25)1 
provide us with 80 hole and particle solutions while only 
16 solutions are obtained with ht=1. Therefore, con- 
tributions to Af(u>) with a more collective nature can be 
better accounted for in the former case (i.e., jit=5). This 
is further corroborated by comparing our DOS, computed 
with riT—5 transformations, with the one obtained using 
an ED, performed with an in-house code, shown in panel 
b) of the figure. Note that we have intentionally used 
a small broadening r=0.05t to retain as much structure 
as possible in our DOS as well as to emphasize the dif- 
ferences with the ED one. As can be observed there is 
excellent agreement in the position and relative heights 
of all the prominent peaks. The hole (blue) and parti- 
cle (black) spectral functions, computed with ny=5 HF- 
transformations, are displayed in panel c) of the same 
figure. We have not included the ones provided by the 
ED since they are quite similar to ours. Their structure 
is dominated by a main peak but less prominent ones 
are also visible in the figure. The momenta (0,7r) and 
(ir, 0), at the noninteracting Fermi surface e(k Q ) = [see, 
Eq. ([A6p of appendix [A] . have the largest spectral weight 
near w-U/2=0. On the other hand, the momenta (0, 0) 
and (0,±7r/2) {(it, ±7t/2) and (tt, tt)] inside (outside) the 
noninteracting Fermi surface contribute mostly to hole 
(particle) states. 
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FIG. 8: Ground state energy of the 4x4 lattice with N e = 15 electrons at U=4t computed with various approaches. The 
different columns refer to the unprojected Hartree-Fock (HF) calculation, HF with linear momentum projection (LM), HF with 
projection of linear momentum and only the z-component of the total spin (LM+S Z ) and HF with projection of linear momentum 
and full spin projection before the variation (LM+S). For all these methods we have used the approximation discussed in Sec. 
Ill Bl with five transformations. Note that the LM+S method corresponds to the symmetry-projected configuration mixing 
approach used throughout the paper. The predicted energies are compared with the exact (EXACT) one.— For details, see the 
main text. 



In Fig. [31 we display the occupation numbers of the ba- 
sis states [see Eq.©] in the 9° = (0, 0, 0) ground state of 
the half- filled 2x4 lattice. Results are shown for the on- 
site repulsions U=4t, 20t, 40t, 64t, 80t and 120t. The cal- 
culations were performed using 80 hole solutions h (i.e., 
71t=5 HF-transformations) in Eq. (l26l) . The evolution of 
the occupations clearly depict the transition to the strong 
coupling regime where the Hubbard Hamiltonian^ can be 
mapped into the AF Heisenberg model In fact, for U 
> 64t the results look very similar to the uniform dis- 
tribution, with occupations = 1, expected in the 
limit U — > oo. 



B. The square 4x4 lattice 

In panel a) of FigJU we show the energies e® obtained, 
via Eq.(gT]), for the half- filled 4x4 lattice at U=4t. Re- 
sults are only shown for the six essentially different pairs 
of linear momentum quantum numbers (0,0), (1,0), (1,1), 
(2,0), (2,1) and (2,2). For each of them we have plotted 
the energies of the first five solutions with spins S=0,1, 
and 2. In this case, the number of variational param- 
eters in our approximation is n var (S — 0,^,^)^512, 
n var (S = l,Cx,C a )=516 and n var (S = 2, £ x , 4)=520 
while the dimensions of the restricted Hilbert spaces are 
n RH (S = O^ty) « 2 x 10 6 , n RH {S = 1, « 
4 x 10 6 and n RH (S = 2,£ B ,£ I/ ) « 3 x 10 6 , respectively. 

The energy -13.5898t of our 6° = (0, 0, 0) ground state 
accounts for 99.76 % of the exact one^^ -13.6219t. In 
order to put our result in perspective, the relative error 
0.24 % in our ground state energy per site ef /16 at U=4t 
should be compared, for example, with the value 0.70 % 



recently reported^ - within the framework of the varia- 
tional MC (VMC) approximation using an ansatz, con- 
sisting of the product of a correlator product state tensor 
network and a Pfaffian wave function, with 524,784 vari- 
ational parameters. Note, that the DMRG formalism in 
momentum spaced (kDMRG) predicts a relative error 
of 0.37%. We have also studied our ground state energy 
per site ef /16 as a function of the interaction strength 
U=2t, 6t, 8t, lOt, 12t and 16t and found relative errors 
always smaller than 0.4 %, 

Coming back to the spectrum shown in panel a) of 
FigHl we note that the first excited state corresponds 
to a 9 = (1,2,2) configuration [with linear momenta 
(tt, 7r)]. In fact, similar to the half-filled 2x4 lattice, 
the first excited state for each combination (£ x ,£y) has 
spin S=l, exception made of (0,0). The 2x4 lattice 
displays a low-lying S=0 singlet (see FigJTJ) while an S=2 
quintet appears in the 4x4 lattice. The excitation en- 
ergies, referred to the 9° = (0, 0, 0) configuration, of 
these low-lying S=l and S=2 states are shown in panel 
b) of Fig|4] as functions of the linear momentum quan- 
tum numbers. The shape of the curve does not fully 
agree with the one obtained with the spin-density wave 
(SDW) approximation 62 mainly due to the absence of de- 
generacy between the r=(0,0) and Q=(2,2) as well as the 
two peaks for the i?i=(l,0) and i?2=(2,l) points. Much 
of this discrepancy could, however, be due to finite size 
effects*^ Note that the two peaks at i?i and i?2, result- 
ing from a kinetic-energy gap of 2t, are already visible 
for the Fermi gas (U=0t). 

The DOS [Eq.pO])] for the half-filled 4x4 lattice 

at U=4t is shown in panel a) of Fig[S] The calculations 
have been performed using ut=5 HF-transformations 
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FIG. 9: (Color online) The energy spectrum, obtained via Eq. (|21[) . for the 4x4 lattice with N e = 14 electrons at U=4t is 
shown in panel a). In panel b), the excitation energies from the ground state to the lowest-lying S=0,1, and S=2 states from 
panel a) are plotted as functions of the linear momentum quantum numbers F = (0,0), Ri = (1,0), P = (2,0), R2 = (2, 1), 
Q = (2,2), and R3 = (1, 1), respectively. In addition to U=4t (blue boxes), results for U=0t (red diamonds) are also included 
for comparison. 



along the lines described in Sec. Ill CI In this case, 
Eqs. (j2~4"| and (|2"51) provide us with 160 hole and particle 
solutions. A Lorentzian folding of width T=0.2t has been 
used. Similar to the case of the half- filled 2x4 lattice, the 
Hubbard gap, Ah = U/2 = 2t, remains present in this 
larger system. The hole (blue) and particle (black) spec- 
tral functions shown in panel b) of the figure show that 
the momenta (±7r/2, ±7r/2), (0, n) and (7r,0) at the non- 
interacting Fermi surface have the largest spectral weight 
near cj-U/2=0. Moreover, the spectral weight due to 
these momenta at the Fermi surface is particle-hole sym- 
metric. 

In panel a) of FigUH we have plotted the DOS Af(ui) 
[Eq.dnO])] for the half-filled 4x4 lattice at U=8t (i.e., 
an on-site repulsion equal to the noninteracting band- 
width W=8t) computed with nx=5 HF-transformations. 
A Lorentzian folding of width r=0.2t has been used. The 
corresponding hole (blue) and particle (black) spectral 
functions are also displayed in panel b) of the figure. 
Our DOS and spectral functions [in particular, the ones 
corresponding to the linear momenta (7r,0), (7t/2, 0) and 
(0, 0)] can be compared with the ones, obtained using the 
Lanczos method, shown in Fig. 2 of Ref. 57. As can be 
observed, the main qualitative features of the particle- 
hole symmetric DOS are well reproduced, namely the 
two prominent peaks at w-U/2 « 2t and 3t (-2t and -3t) 
a lump peaked around oj-U/2 « 5t (-5t) and a smaller 
satellite peak in the neighborhood of w-U/2 sa 8t (-8t). 
In agreement with the results of Ref. 57, the upper and 



lower bands as well as the Hubbard gap are also clearly 
visible in FigO 

We have also studied the evolution of the DOS for the 
half-filled 4x4 lattice as a function of U. To this end, in 
addition to the cases U=4t and 8t shown in Figs. [S]and|ni 
calculations have also been performed for U=2t, 12t and 
20t. In good agreement with previous studies,— ~ 41 ' 57 i 59 
we observe that the Hubbard gap persists for increasing 
values of U. In our calculations, a pronounced suppres- 
sion in the DOS around w-U/2=0 is observed with the 
DOS fully vanishing around U=8t (i.e., around the nonin- 
teracting bandwidth W=8t) which is precisely the region 
where a sizeable Hubbard gap is developedi 42 i 57 ' 59 

Let us now consider two examples of a doped 4x4 lat- 
tice. In FigO we show the spectrum in the case of 15 elec- 
trons at U=4t. For each of the linear momentum quan- 
tum numbers (0,0), (1,0), (1,1), (2,0), (2,1) and (2,2), we 
plot the energies of the first five solutions of Eq. (|2"T]) for 
the spins S=l/2 and 3/2. The number of variational pa- 
rameters in our approximation n var (S = 1/2, £ y )=512 
and n var (S = 3/2, £ x , £ y )=516 should be compared with 
the dimensions tlrh{S = l/2,^ x ,^ y ) 2 xlO 6 and 
n RH (S = 3/2,^,4) « 2 xlO 6 of the restricted Hilbert 
spaces. The first noticeable feature in Fig(7] is that the 
four-fold degenerate O - = (1/2,1,1) ground state has 
non-zero linear momenta (7r/2,7r/2). A finite linear mo- 
mentum for the one-hole ground state has also been pre- 
dicted in previous studies 2 using a variety of approxima- 
tions for lattices of different sizes. Our numerical calcula- 
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FIG. 10: (Color online) The same as Fig[S]but for the 4x4 lattice with iV e = 14 electrons at U=4t. 



tions also predict a two-fold degenerate (1/2,2,0) config- 
uration whose energy is almost the same as the ground 
state one. For the noninteracting system (U=0t), the 
lowest-lying S=l/2 and S=3/2 states with linear momen- 
tum quantum numbers r = (0,0), P = (2,0), Q = (2, 2), 
and i?3 = (1,1) are degenerate and the same is also true 
for the configurations R\ = (1, 0) and R2 = (2, 1). There- 
fore, the huge degeneracy observed in the noninteracting 
case is already partially lifted at U=4t. 

In FigJSJ we compare the ground state energy of the 
4x4 lattice with 15 electrons with the exact oneil 
for U=4t. The energy -14.5469t predicted within our 
symmetry-projected configuration mixing approach, via 
Eq. (l2"Tl) . accounts for 99.19 % of the exact result. It is 
interesting to note that linear momentum plus S z projec- 
tion already accounts for 98.41 % of the exact solution. 
Nevertheless, full spin projection, while also recovering 
the total spin quantum number, still brings a sizeable 
amount of correlations. 

The (shifted) differences ef° - e^ 1/2, ^ :Cy) -U/2, where 
ef is the ground state energy of the half-filled lattice 
(FigfJI and g^ 1 / 2 ^'^ represents the energy of each of 
the lowest-lying S=l/2 states in Figj7l compare very well 
with the position of the first prominent peak in the hole 
spectral functions shown in panel b) of Fig(5] For ex- 
ample, the variational approach predicts ef ° — e^ 2 ' 1 ' 1 ^- 
U/2=-1.044t while the corresponding peak in the hole 
spectral function is predicted to be at w-U/2= -l.OlOt. 
The same is also true for the configuration with linear 
momenta (7r,0) for which the variational approach pre- 
dicts ef° — e^ 2 ' 2 ' ' ) -U/2=-1.052t whereas the position 



of the corresponding peak in the hole spectral function 
is predicted to be at w-U/2= -1.012t. This leads to the 
conclusion that the lowest-lying S=l/2 states in the spec- 
trum of FigfT] are reasonably well described by a wave 
function of the form (|23|) . This is remarkable, since no 
orbital relaxation is accounted for in this wave function, 
i.e., the determinants |Z>W) in Eg. (f2"3")l correspond to the 
ones obtained at half-filling. 

The spectrum obtained, via Eq. (|2"Tj) , for 14 electrons at 
U=4t is displayed in FigJHl The number of variational pa- 
rameters in our approximation n var (S = 0, £, x , £^=504, 
n var (S = 1,6s, £^=508 and n var (S = 2,£ a ,|j / )=512 
should be compared with the dimensions urh{S — 
0,Zz,Sv) « 10 6 , n RH (S = U*,^) PS 2 x 10 6 and 
nRn(S = 2,£ x ,£g / ) « 10 6 of the restricted Hubert spaces. 
The ground state corresponds to the 6°=(0,2,2) config- 
uration [linear momenta (tt, n)] with energy -15.5872t, 
while the exact one is -15.7446tfil On the other hand, 
the VMC approximation 53 predicts a ground state en- 
ergy of -15.5936t. Thus, both methods, ours and VMC, 
yield essentially the same relative error of around 1 % in 
the ground state energy per site. On the other hand, a 
relative error in the ground state energy per site of 0.45 
% is obtained within the kDMRG approximation. 61 

Our calculations also predict two other close-lying 
(0,2,2) solutions (with energies -15.5747t and -15.5777t) 
which cannot be distinguished in Fig|5]and therefore ap- 
pear, together with the ground state, as a single thick 
black line. The energy -15.5743t of the ©=(0,0,0) con- 
figuration is also close to the actual ground state. As a 
result, the symmetry of the T and Q points is almost re- 
covered in panel b) of FigJ5] where the energies of the 
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FIG. 11: (Color online) The same as Fig[5]but for the 4x4 lattice with N e — li electrons at U=8t. The shapes of the DOS as 
well as the spectral functions are qualitatively similar to the ones obtained using Lanczos calculations.— 



lowest-lying S=0,l,2 states for each linear momentum 
combination are shown, referred to the 0° = (0, 2, 2) 
ground state for this system. Note that the configura- 
tions (0,1,1) and (0,2,0) (i.e., the points i? 3 and P in 
panel b) of FiglH]) have very small excitation energies, 
0.0552t and 0.1242t, respectively. The two peaks at the 
points R\ and i?2 are also present in the system with 
N e = 14 electrons at U=0t. The spectrum in panel a) of 
FigjS] exhibits an increase in the density of energy levels, 
compared to the one at half-filling, pointing to its very 
correlated nature. 

The DOS JV(oj) [Eq.flU] for 14 electrons at U=4t is 
shown in panel a) of FigJTU] The calculations have been 
carried out by approximating the (N e ± l)-electron sys- 
tems with ht—5 HF-transformations along the lines de- 
scribed in Sec. Ill CI A Lorentzian folding of width T=0.2t 
has been used. The hole (blue) and particle (black) spec- 
tral functions are displayed in panel b) of the figure. In 
this case, Eqs. (|2"4l and (|2"5)) provide us with 140 hole and 
180 particle solutions. The chemical potential is now lo- 
cated around w-U/2=-1.2t. The comparison with panel 
b) of FigJSJreveals that the structure of the hole states for 
w-U/2 < -2t [i.e., those with linear momenta (±7r/2,0), 
(0, ±7r/2) and (0, 0)] remains to a large extent intact. On 
the other hand, a large fraction of the particle spectral 
weight observed at half-filling for It < cj-U/2 < 3t is 
removed. This depletion occurs in favor of new states 
near w-U/2=0. The spectral decomposition of the DOS 
clearly shows that it is states around the Fermi surface 
that suffer the most pronounced changes with respect to 
half-filling. As a result, the particle-hole symmetry in 



the DOS, observed in FigJSl is suppressed for this doped 
lattice and the original gap dissapears. 

In panel a) of FigfrTl we have plotted the DOS Af(ui) 
[Eq.([2n])] for the 4 x 4 lattice with 14 electrons at 
U=8t. The calculations were performed with 77,^=5 HF- 
transformations and a folding T=0.2t has been used. The 
hole (blue) and particle (black) spectral functions are dis- 
played in panel b) of the figure. Our DOS and spectral 
functions can be compared with the ones, obtained using 
the Lanczos method, shown in Figs. 3 and 4 of Rcf.57. 
The main qualitative features of the DOS are well repro- 
duced, namely the prominent peaks around w-U/2=-4t, 
-3t, -2t and 4t. As can be noted from panel b) of FigfTTl 
the chemical potential is now located around w-U/2=- 
2.4t. One of the main features of the DOS is that it 
displays a pronounced pseudogap which, as can be seen 
from Figs. [6land [TTl results mainly from pulling particle 
strenghts into the (half-filling) gap combined with size- 
able contributions of particle states around w-U/2=4t. 
We note, that the pseudogap problem in the doped 2D 
Hubbard model has also received attention within the 
framework of quantum cluster approaches^ 

We have also performed calculations for the 4x4 lattice 
with 14 electrons at U=2t, 12t and 20t. From these cal- 
culations, and the results already discussed above for the 
cases U=4t and 8t, we conclude that upon dopping with 
two holes the original gap observed at half-filling dissa- 
pears for U=2t and 4t, a pseudogap is developed around 
the noninteracting bandwidth W=8t while for the larger 
interaction strengths U=12t and U=20t the gap is not 
filled. Similar conclusions have been obtained within the 
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FIG. 12: (Color online) The same as FigQ]but for the half-filled 6x6 lattice at U=4t. 

From FigIT2l we realize that the first excited state cor- 
responds to a O — (1,3,3) configuration [with linear 
momenta (7r,7r)]. The energy difference 0.1331t between 
this (it, 7r)-configuration and the ground state is smaller 
than the corresponding value 0.1651t for the half-filled 
4x4 system. On the other hand, similar to the half- 
filled 2x4 and 4x4 lattices, most of the first excited 
states for each combination (£ x , £ y ) have spin S=l, excep- 
tion made of (0, 0) and (2, 3) for which an S=2 quintet 
appears. Note that our calculations predict the lowest- 
lying (1,0,1), (1,1,1) and (1,2,3) solutions to be quite close 
in energy. As with the other half-filled lattices studied, 
we find only a handful of excited states within an energy 
window of t from the ground state. 

The DOS Af(u) [Eq.®] for the half-filled 6x6 lattice 
at U=4t is shown in panel a) of FigfT31 The calculations 
have been carried out with n,T=5 HF-transformations 
along the lines described in Sec. Ill CI In this case, 
Eqs. ([Ml) and (|2"8|) provide us with 360 hole and parti- 
cle solutions. A Lorentzian folding of width T=0.2t has 
been used. The DOS for this large lattice shows a clear 
Hubbard gap Ah = U/2 = 2t. The fact that the gap 
remains intact in going, at half-filling, from the 4 x 4 to 
the 6x6 lattice is consistent with previous studies within 
the DCA which show that it is preserved even in the 
TDLj22t— We note, that the DMFT (i.e., N c = 1) does 
not predict a gap for U=4t. It is only for a larger number 
N c of clusters that the gap starts to develop in these stud- 
ies at the TDL, 37 Back to our DOS, we see again from 



Lanczos framework^ 



C. The square 6x6 lattice 

Finally, let us turn our attention to the half-filled 6x6 
lattice at U=4t. The dimensions ujih{S — 0, £<;,£«) ~ 
2 x 10 17 , n RH (S = l,^,e y ) « 6 x 10 17 and n RH {S = 
2,£x,£j/) ~ 6 x 10 17 of the corresponding restricted 
Hilbcrt spaces are far too large for a brute force diag- 
onalization to be feasible. Other approximate methods 
are then called for, not only to describe ground state 
properties but also to access the excitation spectrum in 
this relatively large lattice for which information is rather 
scarce. In this case, the number of variational parame- 
ters in our approximation is n var (S = 0, £ x > £j/) = 2592, 
n var (S = l,£a!)£tf)=2596 and n var (S = 2,£ !C ,£ 1/ )=2600, 
respectively. Therefore, the half-filled 6x6 lattice repre- 
sents a very challenging testing ground for our symmetry- 
projected configuration mixing approximation. 

In Fig[l2l we show the energies obtained, via 
Ea.(|2ip. for the ten essentially different pairs of linear 
momentum quantum numbers (0,0), (0,1), (0,2), (0,3), 
(1,1), (1,2), (1,3), (2,2), (2,3) and (3,3). For each of them 
we have plotted the first five solutions with spins S=0,1, 
and 2. The ground state corresponds to the 8°=(0,0,0) 
configuration with energy ef =-30.5766t. This can be 
compared with the energy obtained using state-of-the- 
art auxiliary-field MC, -30.89(l)t^ 
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FIG. 13: (Color online) The same as Fig0but for the half-filled 6x6 lattice at U=4t. 



the spectral decomposition, shown in panel b) of FigfTBl 
that it is states at [i.e., (±27r/3,±7r/3), (±7t/3, ±2tt/3), 
(tt, 0) and (0, 7r)] or close to the noninteracting Fermi 
surface that contribute the largest spectral weight to the 
prominent hole and particle peaks around w-U/2=-t and 
w-U/2=t. In general, the DOS for this finite size lattice 
is still highly peaked, though features should smooth out 
as one approaches the TDL. Systematic studies for the 
8x8 and 10 x 10 lattices are in progress and will be 
presented elsewhere. 

IV. CONCLUSIONS 

How to accurately describe many-fermion systems with 
approximate methods, which truncate the complete ex- 
pansion of the wave functions to a numerically feasi- 
ble number of configurations, is a central question in 
nuclear structure theory, quantum chemistry, and con- 
densed matter physics. To this end, in the present study 
we have explored an alternative avenue for the 2D Hub- 
bard model. The main accomplishments of the present 
study are: 

• We have presented a powerful methodology of a 
VAP configuration mixing scheme, originally de- 
vised for the nuclear many-body problem, but not 
yet used to study ground and excited states, with 
well defined quantum numbers, of the 2D Hubbard 
model with nearest-neighbor hopping and PBC. 



Our scheme relies on the Ritz variational principle 
to construct, throughout a chain of VAP calculations, 
a truncated basis consisting of a few (orthonormalized) 
symmetry-projected HF states. The simple structure of 
the projected wave functions employed, combined with a 
fast minimization algorithm, allows to keep low compu- 
tational cost in building our basis. A further diagonal- 
ization of the Hamiltonian within such a basis allows to 
account, in a similar fashion, for residual correlations in 
the ground and excited states. 

• Due to the simple structure of the wave functions 
in our approximation, we can construct an ansatz 
[Eqs. ((221 ) and (|2"?j)]. whose flexibility is well- 
controlled by the number of HF-transformations 
included, to approximate the ground state of the 
(N e ± l)-electron system. This allows us to deter- 
mine one-electron affinities and ionization poten- 
tials as well as to access the spectral weight of states 
with different linear momentum quantum numbers 
in the calculation of spectral functions and the cor- 
responding density of states. 

• We have shown that our approximation gives accu- 
rate results, as compared with exact energies, for 
the 2x4 and 4x4 lattices. We have also provided 
the low-lying spectrum of the 6x6 lattice which, to 
the best of our knowledge, has not been reported 
in the literature. Our ground state energy for this 
lattice compares well with results from state-of-the- 
art auxiliary-field Monte Carlo calculations. 
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Regarding the physics of the 2D Hubbard model, we 
have discussed the trends, in going from the 2 x 4 to 
the 4x4 and 6x6 half-filled lattices, of both the low- 
lying spectra and the spectral functions as well as the 
corresponding density of states. We have found that the 
ground states correspond to configurations with spin zero 
and linear momenta (0,0). We have also found that most 
of the lowest-lying excited states display spin S=l. The 
doped systems with 14 and 15 electrons in the 4x4 lattice 
have also been considered. The ground states of such sys- 
tems correspond to configurations with linear momenta 
different from zero. 

Special attention has been paid to the spectral weight 
of states with different linear momentum quantum num- 
bers. We have compared the DOS predicted within our 
approximation with the one obtained using an exact di- 
agonalization for the half-filled 2x4 lattice and found 
an excellent agreement between the two. Our results for 
the half-filled 4x4 lattice, at different on-site repulsions, 
agree qualitatively well with the ones obtained using the 
Lanczos method. 57 For all the considered half-filled lat- 
tices, a Hubbard gap is predicted within our approxi- 
mation. In particular, the fact that this gap persists in 
going from the 4 x 4 to the 6x6 system is consistent with 
previous studies within cluster extensions to dynamical 
mean field theory which show that it is preserved even in 
the thermodynamic limit. As opposed to the half-filled 
case, the particle-hole symmetry in the DOS is removed 
when doping is present in the system. From the calcu- 
lations for 14 electrons in the 4x4 lattice we conclude 
that for on-site repulsions smaller than the noninteract- 
ing bandwidth the (half-filling) gap dissapears, a pseudo- 
gap develops around U=8t while the gap is not filled for 
larger U values. These results agree well with similar con- 
clusions extracted from Lanczos calculations. 5 - We have 
also found the remarkable result that all the lowest-lying 
S=l/2 states in the spectrum of the 4x4 lattice with 
15 electrons can be reasonably well described by a wave 
function of the form (1231) in which no orbital relaxation 
is accounted for. 

One important feature of the scheme presented in this 
study is that it leaves ample space for further improve- 
ments and research. First, the number of symmetry- 
projected configurations in our basis set can be easily 
increased. Second, we could still incorporate particle 
number symmetry breaking and restoration in our config- 
uration mixing scheme to access even more correlations. 
Our methods could be useful even for more complicated 
lattices like the honeycomb one. An extension of the con- 
sidered 2D t-U Hubbard Hamiltonian to the t—t —t —U 
case is also straightforward, allowing the study of several 
interesting issues like indications of spin-charge separa- 
tion in 2D systems (see, for example, Ref. 65 and refer- 
ences therein). Last but not least, not only the configura- 
tion mixing scheme applied in the present work but also 
the full hierarchy of approximations discussed in Ref. 43 
can be implemented for the molecular Hamiltonian in the 
realm of quantum chemistry, within the already success- 



ful PQT. ■ Work along these avenues is in progress. 

We would like to stress that the cost of the symmetry- 
projected calculations described in this work has the 
same scaling as mean-field methods i 47 ! 48 This statement 
is true as long as the number of grid points required in the 
symmetry restoration remains relatively constant (this is 
usually the case for spin and number projection). Note, 
however, that the restoration of translational symmetry 
in Hubbard lattices with PBC requires a number of grid 
points equal to N S i tes ■ This makes the cost of our calcu- 
lations 0(N S i tes ) more expensive than the Hartree-Fock 
method, which is still a very reasonable scaling. The 
computational effort is mainly concentrated in looping 
over grid points for the evaluation of matrix elements. 
This task is trivially parallelizable and one can thus eas- 
ily reach clusters larger than the ones considered in this 
study. Preliminary calculations for the half-filled and 
doped 8x8 and 10 x 10 lattices are in progress. 

A discussion of the limitations of our method is also in 
order here. Evidently, our method relies on the Hamil- 
tonian having good symmetries. The lower the number 
of symmetries of a given Hamiltonian, the lower the cor- 
relations that can be accounted for by means of symme- 
try restoration. On the other hand, the most interesting 
quantum behavior is found in systems where symmetries 
are present. Second, it is known that the correlation 
energy per particle obtained with approaches based on 
a single symmetry-projected determinan t 47 ! 48 decays as 
the lattice size increases. We observe that the error in 
the energy per site is larger in the case of the half-filled 
6x6 lattice than in the 4x4 one, although in both cases 
we have restricted the present study to m—5 transforma- 
tions. One can, however, increase the number of trans- 
formations to maintain the quality of our wave functions. 
In principle, if the number of transformations is equal to 
the size of the restricted Hilbert subspace the method 
becomes exact. In practice, one can only hope that the 
number of transformations needed to access the relevant 
physics of the considered lattices is relatively low. 

Finally, we believe that the finite size calculations dis- 
cussed in the present work are complementary to other 
approaches where impurity solvers play an important 
role. However, the symmetries to be broken and restored 
in the impurity-bath and bath Hamiltonians will depend 
on the details of the case considered. 42 
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Appendix A: Symmetry-projected matrix elements between two Slater determinants \D l ) and \T> k ) 

In this appendix, we present the expressions for the matrix elements 70£® = (T> l \HH U i ) P® T .,\'D k ) and jV^w = 
/V l \P®, \ v k\ requ i rec i to com pute the kernels H™3 and Af^S in Eq.® of Sec. iHBl Note that the matrix elements 

2-i2-i ____ 2-i2-i ZjZj 

required in Eq.® of Sec. Ill Al are just a particular case where both Slater determinants are the same. Here, and in 
what follows, we keep our notation as close as possible to the one already used for the ID Hubbard model. 44 Both 
and Af*3 read 



ft' V V 7 



ES ' " 8ir 2 N 



sites 

J 



rike 25 + 1 



jvnfcfc) 



J 



^ e"^ J cMD^, (Q)h ik (Q,j)n ik (Q,j) 
j 

]T J dflD^, (Q)n ifc (fi, j) (Al) 



where, k^ = (fc^ x , fc^) = (t^ 2 -, "75^") anc ^ J = (jx,jy), respectively. For the gauge-rotated norm 

n ik {Sl,j)=det Ne X ik (n,j) (A2) 
the determinant has to be taken over the N e x N e dimensional occupied part of the matrix 

X$(Sl t j) = {-D* T S(n,j)V k *) ab (A3) 

with 

s aaa ,(n,i)=v 1 J*(ny^ (A4) 

The gauge-rotated Hamiltonian takes the form 

v k (n,j) = if fe (o, j) + ii> (r ifc (nj)p w (o, j)) (A5) 

with 

e(k Q ) = -2t (cosk ax + cosfe ay ) 



ly sitr- 



X 



/35 

x '5 aa>P t_^ p _^)-{l-Ka')pt,p-A^S)\ (A6) 

where the product of generalized Kronecker deltas in T lk (n,, j) results from the transformation of the on-site interaction 
term in Eq.([T]) to the momentum representation. As a consequence of the PBC, S^^^. _ s , is one if on + pi — 7j — <5j 
is either or ±JVj and zero else. 

Appendix B: Symmetry-projected particle-hole matrix elements between two Slater determinants |2? 1 ) and 

\V k ) 

In this appendix, we present the expressions for the matrix elements % tk ®,' ph = l(D % \HHubP® Y: i bj,('D k )bh('D k )\'D k ) and 
tfik®-, P h = (23*|^ E ,St(2?*)Sfc(23 fc )|23*> required to compute the kernels K,™°; ph and 1l™®; ph defining the variational 
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equations discussed in Sec. Ill Bl Note that the matrix elements required in Eq. (flQ|) of Sec. Ill Al are just a particular 
case where both Slater determinants arc the same. We obtain 



ik®\ph 

BE' ~~ 



2S+1 
25+1 



£ e-^J J dQD 3 ^, (OK fe (fi, (£1, j) 



= 8^^E e ^ j /^BB'(^K fe («»j)<.(^j) 



(Bl) 



where, as in appendix^ = (fc^ x ,fc^ H ) = f^j^S ~jf\ and j = (j x ,jy), respectively. On the other hand, 



n ik 
n ph 



(n.J)= E [^*'(n.J)] ^(nj) 



(B2) 



with the indices /i (p) running over all the occupied (unoccupied) states in \T> k ). The inverse <¥**,(0,j) 
over the occupied part of the matrix (|A3|) . Finally, 



is taken 



hp 



(B3) 



with the functions y[Ct,j) and W (fi>j)) defined, for all the occupied h and unoccupied p states in \V k ), as 



2?' 



ft' 

Vv, p (nj) = E [i-p^'Cn.j) 



/ ~ „ Oct a V ' J / So- ,p 



<5er" er'" 



"yer ,5cr 



(B4) 



Appendix C: Symmetry-projected matrix elements between two Slater determinants \D l ) and \T> k ) for 

spectral functions 

In this appendix, we present the computation of the kernels W and Af required in Eq. (l24l) as well as of the 
kernels H e+ and Af e+ in Eq.flU). Both W e ~ and A/" e ~ read 



TV" 6 " , , 

ihaikh c 



ihcrikh o 



8tt 2 N s 



j 

2 e- <k ^ j J dfiD^*(n)n tt (n, , (fi, j) 



with the vector = ^fc^-,fc^-^ = f 2 ^ x , j while i,k = 1, ... utj ^, ^ = l,-.-A^ e and cr,<T 



other hand 



(CI) 

±1/2. On the 



°hh' 



(oj)= %(nj 



and 



(C2) 



(C3) 
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respectively. The function h (£l,j) reads 



while y*? a<T (to,i) is given in Eq.(|B4|. 

The norm and Hamiltonian overlaps in Eq. ([28[) read 



(C4) 



a/-: 



ip<7,kp a jy_ 



i£es . J 
J 

ites . J 



with the vector = ^fc^+,fc^+^ = ^ 
On the other hand, 



) while i,k = 1, .. . n T , p,p = N e + 1, . . . 2N sltes and cr, cr' 



(C5) 
±1/2. 



aft' 



(C6) 



and 



*&(nj) = n«,(nj)/ l «(n,j)+ vr(n,j)r*(nj)viHnj) 



life / 



(C7) 



respectively. The function Wp^ aCT (0, j) is given by 



M&a(nj) = E 2 V J ,[ 1 -^(°.J) 



/3<r 



/3cT ,0£(T 



(C8) 



-A'/ 



while W 7<T ' )P '(n,j) is defined in Eq.(|B4 
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